% The economic effects of international administrations: The cases of Kosovo and East Timor
% by C�sar Urquizo and Diego Winkelried (Universidad del Pac�fico, www.up.edu.pe)
% Economic Development and Cultural Change (2019)
% https://doi.org/10.7910/DVN/ENX8VL

clear all
clc

%% Load data
load DataUW2019

[country_list, j] = unique(country);
cnames = Cnames(j);
nc = length(country_list);  % Number of countries
T  = size(Ydata,1)/nc;      % Number of periods per country

% Treated and excluded from donor pool
treated  = 'Kosovo';
excluded = {'Timor-Leste' 'Iraq' 'Kosovo' 'Liberia' 'Eritrea'};
Tdate = 2000;                         % Treatment date
T0 = Tdate-1990+1; 
dates = [T0-10 T0-5 T0];              % Data starts in 1990
excludefe = {'life_exp' 'fert_rate'}; % Exclude covariates for which Kosovo is an outlier

% Pretreatment averages ("fixed effects")
k  = length(Xfenames);
Xfe  = nan(k, nc);
presample = (T0 - 10):(T0 - 1);
for i = 1:nc
   ind = (country == i);
   Xfei = Xdata(ind, :);
   Xfe(:, i) = nanmean( Xfei(presample,:) );
end

% Filtering of outcome data
logY = 100*log(Ydata);
Yhp = nan(size(Ydata, 1), 1);
for i = 1:nc
    ind = (country == i);
    [Yhp(ind), ~] = hpfilter(logY(ind), 6.25);
end

% Pretreatment outcome values
nd = length(dates);
Ypre = nan(nd, nc);
Yname = 'ln_gni_pcft';
for i = 1:nd
    Yprenames{i,1} = [Yname '(' num2str(1989 + dates(i)) ')']; %#ok<SAGROW>
end
for i = 1:nc
    ind = (country == i);
    Yi = Yhp(ind, :);
    Ypre(:, i) = Yi( dates );   
end

% Selecting covariates
if ~isempty(excludefe)
    ok = false;
    drop = false(k, 1);
    for j = 1:length(excludefe)
        i = find(strcmp(Xfenames, excludefe{j}));
        if ~isempty(i)
            drop(i) = true;
            ok = true;
        end
    end
    if ok
        Xfenames(drop) = [];
        Xfe(drop, :) = [];
    end
end

% Selecting the donor pool
controls = cnames;
for i = 1:length(excluded)
    ind = strcmp(controls, excluded{i});
    controls(ind) = [];
end

% Final data
n = length(controls);
Y0 = nan(T, n);
X0 = nan(size(Xfe,1) + size(Ypre,1), n);
Xnames = [Xfenames; Yprenames];
for i = 1:length(controls)
    ind_large = strcmp(Cnames, controls{i});
    Y0(:, i) = Yhp(ind_large);
    
    ind_small = strcmp(cnames, controls{i});
    X0(:, i)  = [Xfe(:, ind_small); Ypre(:, ind_small)];
end
ind_large = strcmp(Cnames, treated);
Y1        = Yhp(ind_large);    
ind_small = strcmp(cnames, treated);
X1 = [Xfe(:, ind_small); Ypre(:, ind_small)];

% For scaling graphs
Ybase = Y1(T0) - 100;

Xnom = char(Xnames);

%% Pretreatment sample
pts = 1:10;

T = size(Y1, 1); 
[K, n] = size(X0);
Kfe = K - 3; % Fixed-effect predictors 


%% Matlab computations
% We optimize for w, but take V as given (this comes from the Stata run)
filename = 'matrix_V_kosovo.txt';
v = loadV(filename, Xnames);
R1 = SyntheticControl(Y1, Y0, X1, X0, pts, v);

%% Sensitivity analysis excluding one fixed-effect predictor at a time
% We renormalize the V matrix
if 01
    Rbase = R1; % Base run
    v = Rbase.v;
    K = length(Rbase.Xs);
    
    YsK = nan(T, Kfe);
    wK  = nan(n, Kfe);
    VK  = zeros(K, Kfe);
    MSEK = nan(Kfe,1);
    dK   = nan(Kfe,1);
    XsK = nan(K,Kfe);
    
    %Loop excluding one predictor at a time
    for i = 1:Kfe
        indK = true(K, 1);
        indK(i) = false;
        
        temp = SyntheticControl(Y1, Y0, X1(indK, :), X0(indK, :), pts, v(indK));
        YsK(:, i) = temp.Ys;
        wK(:, i)  = temp.w;
        MSEK(i)   = temp.MSE;
        XsK(:,i)  = X0*temp.w;
    end
    
    % Catchy graph
    years = 1989 + (1:T)';
    figure

    
    xlim([min(years)-1/2 max(years)+1/2])
    hold off
    
    indK = true(K, 1);
    
    % Upon exploration some X are removed 
    indK([4 5 10]) = false;
    
    Rw = SyntheticControl(Y1, Y0, X1, X0, pts, v, mean(wK,2));
   
    % Report
    BigW   = [Rbase.w  wK Rw.w];
    BigXs  = X0*BigW;   
    BigMSE = [Rbase.MSE MSEK(:)' Rw.MSE];
    BigYs  = Y0*BigW;

    % Data normalization (cross-sectionally)
    WW = zscore([X0 X1]')';
    X0norm = WW(:, 1:n);
    X1norm = WW(:, end);
    Xsnorm = X0norm*BigW;   
    
    %Tables with results
    ncols  = size(BigXs, 2);
    XMnom = strvcat(Xnom, 'Pretreatment MSE', 'Distance of X');

    titl = 'Means of Treated / Synthetic controls and Pretreatment MSE (normalized data)';
    fprintf(2,'\n%s\n%s\n', titl, repmat('-', 1, length(titl)))
    for i = 1:K
        fprintf(1, ['%s \t %8.2f | ' repmat('\t%8.2f', 1, ncols) '\n'], XMnom(i,:), X1norm(i), Xsnorm(i,:));
    end

    titl = 'Weights in synthethic control';
    fprintf(2,'\n%s\n%s\n', titl, repmat('-', 1, length(titl)))
    Cnames = char(controls);
    indsyn = find(sum(BigW, 2) > 1e-5);
    ncols  = size(BigW, 2);
    for i = 1:length(indsyn)
        j = indsyn(i);
        fprintf(1, ['%s' repmat('\t%4.1f', 1, ncols) '\n'], Cnames(j,:), 100*BigW(j,:));
    end
    fprintf(1, '\n');
    
    % Sentitivity analysis graph
    years = 1989 + (1:T)';
    figure
    hold on
    h = plot(years, Y1 - Ybase); set(h(1), 'Linestyle', '-', 'Color', 'b', 'Linewidth', 1); 
    h = plot(years, YsK  - Ybase); set(h(:), 'Linestyle', ':', 'Color', 'r', 'Linewidth', 1);
    [~,j] = min(MSEK);
    set(h(j), 'Linestyle', '--', 'Color', 'r', 'Linewidth', 2);
    h = plot(years, Rbase.Ys - Ybase); set(h(1), 'Linestyle', '-', 'Color', 'r', 'Linewidth', 1);
    h = plot(years, Rw.Ys - Ybase); set(h(1), 'Linestyle', '-', 'Color', 'g', 'Linewidth', 2);
    
    title('\bfVarious runs')

    xlim([min(years)-1/2 max(years)+1/2])
    hold off
end



